clear all
*cap log close
set more off
set seed 603



// ---------------------------------------------------------------------------------------------
// extrapolation for 30+ speeding ------------
local y = "speed_30"
local Q = 2

// store estimates (poly) -----
import delimited using "${est}/speeds_q2.csv", clear
keep if yvar == "`y'"
replace par = trim(upper(par))
keep if par == "Y0"|par=="Y1"
gen Z = 0
replace Z = 1 if par=="Y1"

keep Z est ub lb  
gen flag_poly = 1
tempfile est1
save    `est1'

// store estimates (bins) --------
import delimited using "${est}/speeds_b5.csv", clear
keep if yvar == "`y'"
replace par = trim(upper(par))
keep if par == "Y0"|par=="Y1"
gen Z = . 
replace Z = 0.025 if par=="Y0"
replace Z = 0.975 if par=="Y1"
keep Z est ub lb  
gen flag_bin = 1
tempfile est2
save    `est2'

// estimate quadratic fit ---------
use "${temp}/speeds_data", clear 
forval q = 1/`Q' {
	gen Z_`q' = Z^`q'
}
reghdfe `y' Z_*, absorb(totfe)
local b0 = _b[_cons]
forval q = 1/`Q' {
	local b`q' = _b[Z_`q']
}

// estimate binscatter ---------
use "${temp}/speeds_data", clear 
binsreg `y' Z, absorb(totfe) binspos(es) nbins(12) savedata(bs) replace 

// compile together ------------
use bs, clear 
rename dots_x Z 
rename dots_fit est 
gen flag_bs = 1 
rm bs.dta
qui append using `est1'
qui append using `est2'

// generate quadratic fit for graphing ---------
gen fit = `b0'
forval q = 1/`Q' {
	replace fit = fit + `b`q''*Z^`q'
}

// build plot ----------------
sort Z
#delimit ;
scatter est Z if flag_bs==1, msymbol(Oh) mcolor(dknavy) msize(medium) || 
line fit Z, lcolor(dkgreen) ||
scatter est Z if flag_poly==1, mcolor(dkgreen) msize(med) ||
rspike ub lb Z if flag_poly==1, lcolor(dkgreen) ||
scatter est Z if flag_bin==1, msymbol(D) msize(medium) mcolor(cranberry)
graphregion(color(white)) plotregion(lcolor(black) lwidth(medthin)) 
ylab(,nogrid) xlab(,nogrid)
xtitle("Officer Stringency") ytitle("")
legend(ring(0) cols(1) pos(2) order(1 2 5) region(lstyle(border))
lab(1 "Binscatter") lab(2 "Quadratic Fit") lab(5 "Tail Mean: 5%")) ;
#delimit cr
graph export "${out}/apx_extrap/speed30_extrap.pdf", replace 
// ---------------------------------------------------------------------------------------------





// ---------------------------------------------------------------------------------------------
// levels and gains plot ----------------------------
import delimited using "${est}/speeds_q2.csv", clear
keep if yvar == "`y'"

/*
// For printing ------
qui summ est if par == "Y0D1" 
local b = r(mean)
qui summ se  if par == "Y0D1"
local s = r(mean)
local pt1 = "E(Y0) = `:di %4.3f `b'' (`:di %4.3f `s'')"
qui summ est if par == "ATT" 
local b = r(mean)
qui summ se  if par == "ATT"
local s = r(mean)
local pt2 = "ATT = `:di %4.3f `b'' (`:di %4.3f `s'')"
// For printing ------
qui summ est if par == "Y0D0" 
local b = r(mean)
qui summ se  if par == "Y0D0"
local s = r(mean)
local pu1 = "E(Y0) = `:di %4.3f `b'' (`:di %4.3f `s'')"
qui summ est if par == "ATU" 
local b = r(mean)
qui summ se  if par == "ATU"
local s = r(mean)
local pu2 = "ATU = `:di %4.3f `b'' (`:di %4.3f `s'')"
*/

* populate levels and gains ----------------
do "${code}/utility/create_levels_gains.do"


drop if mi(y0)
replace ind = 0.2 if _n == 1
replace ind = 0.5  if _n == 2 
replace ind = 0.8 if _n == 3
#delimit ; 
twoway 
(rbar y1u y1l ind, color(maroon%20) barwidth(0.02))
(scatter y0 ind, msize(large) msymbol(Oh) mcolor(navy))
(pcarrow y0 ind y1 ind, lcolor(maroon) lwidth(medthick) mcolor(maroon) barbsize(medlarge)),
xscale(r(0 1))
xlabel(0.2 "All Motorists" 0.5 "Treated" 0.8 "Untreated", nogrid) ylab(,nogrid) 
legend(ring(0) pos(2) order(2 3) cols(1) region(lstyle(border))
label(2 "Untreated Level") label(3 "Treatment Effect")) 
xtitle("") ytitle("Pr(Reoffend)") 
graphregion(color(white)) bgcolor(white) plotregion(lcolor(black) lwidth(medthin)) ;
*text(0.31 0.2 "Treated:" "`pt1'" "`pt2'", place(n) size(medium))
*text(0.29 0.2 "Unreated:" "`pu1'" "`pu2'", place(n) size(medium));
#delimit cr
graph export "${out}/apx_extrap/speed30_est.pdf", replace
// ---------------------------------------------------------------------------------------------





// ---------------------------------------------------------------------------------------------
// other plots: by speed 

import delimited using "${est}/speeds_q2", clear 
split yvar, parse("_") gen(stub_)
destring stub_2, gen(speed) force 

gen temp1 = est if par == "ATT"
bysort speed: egen temp2 = mean(temp1)
gen est_scale = est/temp2

#delimit ;
scatter est speed if par=="ATT", msymbol(O) mcolor(midblue%60) yaxis(2) || 
scatter est speed if par=="ATU", msymbol(S) mcolor(maroon%60) yaxis(2)  ||
line est_scale speed if par=="ATU", yaxis(1) lcolor(black) lpattern(dash)
graphregion(color(white)) bgcolor(white) plotregion(lcolor(black) lwidth(medthin))
legend(ring(0) pos(4) order(1 2 3) cols(1) region(lstyle(border))
label(1 "ATT") label(2 "ATU") lab(3 "ATU/ATT (Right Axis)")) 
xtitle("") ytitle("", axis(1)) ytitle("", axis(2))
xlab(,nogrid) ylab(,nogrid) xtitle("Reoffending Speed") ;
#delimit cr 
graph export "${out}/apx_extrap/speeds_gains.pdf", replace

drop temp1 temp2 est_scale 
gen temp1 = est if par == "Y0D0"
bysort speed: egen temp2 = mean(temp1)
gen est_scale = est/temp2

#delimit ;
scatter est speed if par=="Y0D0", msymbol(O) mcolor(midblue%60) yaxis(2) || 
scatter est speed if par=="Y0D1", msymbol(S) mcolor(maroon%60) yaxis(2)  ||
line est_scale speed if par=="Y0D1", yaxis(1) lcolor(black) lpattern(dash)
graphregion(color(white)) bgcolor(white) plotregion(lcolor(black) lwidth(medthin))
legend(ring(0) pos(12) order(1 2 3) cols(1) region(lstyle(border))
label(1 "E(Y0|D=0)") label(2 "E(Y0|D=1)") lab(3 "E(Y0|D=1)/E(Y0|D=0) (Right Axis)")) 
xtitle("") ytitle("", axis(1)) ytitle("", axis(2))
xlab(,nogrid) ylab(,nogrid) xtitle("Reoffending Speed") ;
#delimit cr 
graph export "${out}/apx_extrap/speeds_levels.pdf", replace
// ---------------------------------------------------------------------------------------------


// remove temp data -----
* "rm ${temp}/speeds_data.dta"






